Fast algorithms for computing isogenies 
between elliptic curves * 



Alin Bostan, Bruno Salvy, Projet ALGO, INRIA Rocquencourt 
Domaine de Voluceau, 78153 Le Chesnay Cedex, FRANCE 
{ Alin. Bostan, Bruno. Salvy } Oinria.fr 

Frangois Morain, Eric Schost, LIX, Ecole polytechnique 
91128 Palaiseau, France 
{ Francois .Morain, Eric. Schost } (§lix. polytechnique.fr 

February 1, 2008 



Abstract 

We survey algorithms for computing isogenies between elliptic curves defined 
over a field of characteristic either or a large prime. We introduce a new algorithm 
that computes an isogeny of degree i {I different from the characteristic) in time 
quasi-linear with respect to i. This is based in particular on fast algorithms for 
power series expansion of the Weierstrass p- function and related functions. 

1 Introduction 

In the Schoof-Elkies-Atkin algorithm (SEA) that computes the cardinality of an elliptic 
curve over a finite field, isogenies between elliptic curves are used in a crucial way (see 
for instance jS] and the references we give later on). Isogenies have also been used to 
compute the ring of endomorphisms of a curve [HI] and isogenies of small degrees play a 
role in |23 E] • More generally, in various contexts, their computation becomes a basic 
primitive in cryptology (see |^ H Ell 123 EHl E31 El)- 

An important building block in Elkies's work is an algorithm that computes curves 
that are isogenous to a given curve E. This block uses modular polynomials to get the list 
of isogenous curves and Vein's formulas to get the explicit form of the isogeny I : E ^ E, 
where ii^ is in a suitable Weierstrass form. 

In this work, we concentrate on algorithms that build the degree i isogeny / from 
E and E (and possibly some other parameters, see below). For the special case i = 2, 
formulas exist jlOj; see also ^B]- We could restrict further to the case when i is an odd 
prime, since isogenies can be written as compositions of isogenies of prime degree, the 

*This work was supported in part by the French National Agency for Research (ANR Gecko). 
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case of prime powers using isogeny cycles [THl CEl ESj • Besides, the odd prime case is the 
most important one in SEA. However, our results stand for arbitrary £. 

We demand that the characteristic p of the base field K be or p ^ £. This restriction 
is satisfied in the case of interest in the application to the SEA algorithm, since otherwise 
p-adic methods are much faster and easier to use Several approaches to isogeny 

computation are available in small characteristic: we refer to ^SlISH] for an approach via 
formal groups, [SEj for the special case p = 2, and jTH UHl EH IHI] for the general case of 
p small. The case of p = i deserves a special treatment, see [HI EH], using Gunji's work 
P7j as main ingredient (see also P7j). 

Our assumption on p implies that the equations of our curves can be written in the 
Weierstrass form 

= x^ + Ax + B. (1) 

In characteristic zero, the curve dH) can be parameterized by {x^y) = {p{z) , p' {z) / 2) in 
view of the classical differential equation 

p'izr = Aipizr + Apiz)+B) (2) 

satisfied by the Weierstrass p-function. This is the basis for our computation of isogenies. 
We thus prove two results, first on the computation of the Weierstrass p-function, and 
then on the computation of the isogeny itself. 

Our main contribution is to exploit classical fast algorithms for power series com- 
putations and show how they apply to the computation of isogenies. We denote by 
M : N — > N a function such that polynomials of degree less than n can be multiplied 
in M(n) base field operations. Using the fast Fourier transform ^1 one can take 
M(n) G O (n log n log log n) ; over fields containing primitive roots of unity, one can take 
M(n) e O(nlogn). We make the standard super-linearity assumptions on the function 
M, see the following section. 

Theorem 1. Let be a field of characteristic zero. Given A and B in K, the first n 
coefficients of the Laurent expansion at the origin of the function p defined by ^ can be 
computed in 0{M{n)) operations in K. 

In ^ we give a more precise version of this statement, that handles the case of fields 
of positive, but large enough, characteristic. 

An isogeny is a regular map between two elliptic curves that is also a group morphism. 
If E and E are in Weierstrass form and I = {Ix,Iy) is an isogeny E ^ E, then Ix{P) 
depends only on the x-coordinate of P, and there exists a constant c G K such that 
ly = cyl'^. Following Elkies 1221; consider only so-called normalized isogenies, 
those for which c = 1 (such isogenies are used for instance in SEA). In this case, we will 
write a for the sum of the abscissas of non-zero points in the kernel of /. 

Theorem 2. Let be a field of characteristic p and let E and E be two curves in 
Weierstrass form, such that there exists a normalized isogeny I : E ^ E of degree i. 
Then, one can compute the isogeny I 

1. in 0(M(£)) operations in K, if p = or p> 2^ — 1, if a is known; 

2. in 0(M(£) \ogt) operations in K, if p = or p > 8£ — 5, without prior knowledge 
of a. 
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Taking M(?7,) G 0(?7. lognloglogn) shows that the complexity results in Theorems Q 
and El are nearly optimal, up to polylogarithmic factors. Notice that the algorithms using 
modular equations to detect isogenics yield the value of a as a by-product. However, in 
a cryptographic context, this may not be the case anymore; this is why we distinguish 
the two cases in Theorem |21 

This article is organized as follows. In §21 we recall known results on the fast compu- 
tation of truncated power series, using notably Newton's iteration. In §3 "we show how 
these algorithms apply to the computation of the p-function. Then in ^ we recall the 
definition of isogenics and the properties we need, and give our quasi-linear algorithms; 
examples are given in ^ In the next section, we survey previous algorithms for the com- 
putation of isogenics. Their complexity has not been discussed before; we analyze them 
when combined with fast power series expansions so that a comparison can be made. 
Finally, in ^ we report on our implementation. 



2 A review of fast algorithms for power series 

The algorithms presented in this section are well-known; they reduce several problems 
for power series (reciprocal, exponentiation, . . . ) to polynomial multiplication. 

Our main tool to devise fast algorithms is Newton's iteration; it underlies the 0{M{i)) 
result reported in Theorem d and in the (practically important) point (1) of Theorem |21 
Hence, this question receives most of our attention below, with detailed pseudo-code. We 
will be more sketchy on some other algorithms, such as rational function reconstruction, 
referring to the relevant literature. 

We suppose that the multiplication time function M is super-linear, i.e., it satisfies 
the following inequality (see, e.g., [23 Chapter 8]): 

M(n) M(n') , , 

< if n < n'. 3 

n n' 

In particular. Equation Q implies the inequahty 

M(l) + M(2) + M(4) + ■ ■ ■ + M(2^) < 2M(2^), 

which is the key to show that all algorithms based on Newton's iteration have com- 
plexity in 0(M(n)). Cantor and Kaltofen [10 have shown that one can take M(n) in 
0(?T, lognloglogn); as a byproduct, most questions addressed below admit similar quasi- 
linear estimates. 



2.1 Reciprocal 

Let / = E.>o fi^' be in K[[z]], with /o ^ 0, and let g = l/f = E^>o9i^' K[[^]]. The 
coefficients gi can be computed iteratively by the formula 

1 1 * 

go = — and gi = -— fjgi-j for i > 1. 
Jo Jo ~[ 

For a general /, the cost of computing 1// mod z"' with this method is in O(n^); observe 
nevertheless that if / is a polynomial of degree d, the cost reduces to 0{nd). 
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To speed up the computation in the general case, we use Newton's iteration. For 
reciprocal computation, it amounts to computing a sequence of truncated power series hi 
as follows: ^ 

ho = — and /ij+i = hi{2 — fhi) mod 2;^'^^ for z > 0. 
Jo 

Then, = l/fmodz^\ As a consequence, l//modz" can be computed in 0{M{n)) 
operations. This result is due to Cook for an analogous problem of integer inversion |12j, 
and to Sieveking |1H1 and Kung jHEl in the power series case. 

2.2 Exponentiation 

Let / be in K[[z]], with /(O) = 0. Given n in N, such that 2, ... ,n — 1 are units in K, 
the truncated exponential exp„(/) is defined as 

n— 1 

exp„(/) = J2 if' 
i=o ^■ 

Conversely, if 5? is in 1 + 2;K[[2;]], its truncated logarithm is defined as 

n— 1 

logJ^7) = -5^-(l-^7rmod z\ 

i=l 

The truncated logarithm is obtained by computing the Taylor expansion of g'/g modulo 
z^~^ using the algorithm of the previous subsection, and taking its antiderivative; hence, 
it can be computed in 0{M{n)) operations. 

Building on this. Brent 6j introduced the Newton iteration 

go = 1, gi+i = gi{l + f - logs^+i (s-i)) mod z"^^^' 

to compute the sequence gi = exp2i(/). As a consequence, exp^(/) can be computed in 
0{M{n)) operations as well, whereas the naive algorithm has cost O(n^). 

As an application, Schonhage [IHI gave a fast algorithm to recover a polynomial / of 
degree n from its first n power sums Pi, ■ ■ ■ ,Pn- Schonhage's algorithm is based on the 
fact that the logarithmic derivative of / at infinity is the generating series of its power 
sums, that is. 

Hence, given pi,...,p„, the coefficients of / can be recovered in time 0{M{n)). This 
algorithm requires that 2, . . . , n be units in K. 

2.3 First-order linear differential equations 

Let a, b, c be in K[[2]], with a(0) 7^ 0, and let a be in K. We want to compute the first n 
terms of / G K[[2;]] such that 

(^f + bf = c and /(O) = a 
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Let B = b/a mod ^ and C = c/a mod ^ Then, defining J = exp„(/ B), f satisfies 
the relation 

f = i(^a + j CJ^ mod z". 

Using the previous reciprocal and exponentiation algorithms, / mod can thus be com- 
puted in time 0(M(n)). This algorithm is due to Brent and Kung [Hj; it requires that 
2, . . . , n — 1 be units in K. 

2.4 First-order nonlinear differential equations 

We only treat this question in a special case, following again Brent and Kung's article |Hl 
Theorem 5.1]. Let G be in K[[z]][t], let a, (3 be in K, and let / G K[[2;]] be a solution of 
the equation 

f^ = G{z,f), fiO) = a, f'i0)=/3, 

with furthermore = G{0, a) ^ 0. Supposing that, for s > 2, the initial segment /i = 
/ mod is known, we show how to deduce / mod z"^^'^. Write / = /i + /2 mod z"^^'^, 
where z'^ divides f2- One checks that /2 is a solution of the linearized equation 

2f[f^ - Gt{z, h)h = G{z, A) - /f mod z'^-\ (4) 

with the initial condition /2(0) = 0, where Gt denotes the derivative of G with respect to 
t. The condition /'(O) ^ implies that f[ is a unit in K[[2]]; then, the cost of computing 
/2 mod z^'^"^ is in 0(M(s)) (remark that we do not take the degree of G into account). 
Finally, the computation of / at precision n is as follows: 

1. Let f = a + jSz mod z"^ and s = 2; 

2. while s < n do 

(a) Compute / mod z^^~^ from / mod z'^] 

(b) Let s = 2s-l. 

Due to the super- linearity of M, / mod 2" can thus be computed using 0{M{n)) opera- 
tions. Again, we have to assume that 2, . . . , n — 1 are units in K. 

2.5 Other algorithms. 

We conclude this section by pointing out other algorithms that are used below. 

Power series composition. Over a general field K, there is no known algorithm of 
quasi-linear complexity for computing f{g) mod z", for f,g in K[[2;]]. The best results 
known today are due to Brent and Kung [Hj. Two algorithms are proposed in that article, 
of respective complexities 0{M{n)^Jn + and O ( M (n) y/n log n) , where 2 < < 3 is 

the exponent of matrix multiplication (see, e.g., Chapter 12]). Over fields of positive 
characteristic p, Bernstein's algorithm for composition "3] has complexity 0(M(n)), but 
the 0{ ) estimate hides a linear dependence in p, making it inefficient in our setting 
(p > n) . 



5 



Rational function reconstruction. Our last subroutine consists in reconstructing 
a rational function from its Taylor expansion at the origin. Suppose that / is in K(2) 
with numerator and denominator of degree bounded respectively by n and n', and with 
denominator non-vanishing at the origin; then, knowing the first n + n' + 1 terms of the 
expansion of / at the origin, the rational function / can be reconstructed in 0(M(n + 
n') log(n + n')) operations, see [Zj. 



3 Computing the Weierstrass p-function 
3.1 The Weierstrass p-function 

We now study the complexity of computing the Laurent series expansion of the Weier- 
strass p-function at the origin, thus proving Theorem ^ We suppose for a start that the 
base field K equals C; the positive characteristic case is discussed below. Let thus A, B 
be in K = C. The Weierstrass function p associated to A and 5 is a solution of the 
non- linear differential equation its Laurent expansion at the origin has the form 

p{z) = ^^+Y,c,z^\ (5) 

i>l 

The goal of this section is to study the complexity of computing the first terms Ci, . . . , c„. 
We first present a "classical" algorithm, and then show how to apply the fast algorithms 
for power series of the previous section. 



3.2 Quadratic algorithm 

First, we recall the direct algorithm. Substituting the expansion into Equation 
and identifying coefficients of z'"^ and gives 

ci = -— and C2 = - — . 

5 7 

Next, differentiating Equation ((2|) yields the second order equation 

p" = + 2A. (6) 
This equation implies that for /c > 3, is given by 

3 ^'"^ 

" (fc-2)(2A; + 3) ^ ^^^^-1-- ^'^) 

Hence, the coefficients ci, . . . , can be computed using 0{n^) operations in K. 

If the characteristic p of K is positive, the definition of p as a Laurent series fails, 
due to divisions by zero. However, assuming p > 2n + 3, it is still possible to define the 
coefficients Ci, . . . , c„ through the previous recurrence relation. Then, again, ci, . . . , c„ 
can be computed using 0{n'^) operations in K. 
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3.3 Fast algorithm 

We first introduce new quantities, that are used again in the next section. Define 
Q{z) = ^ e + z^'Kiiz^]] and Riz) = ^/Qi^ E z + z'K[[z% 

The differential equation satisfied by R is 

R'{z)^ = BR{zf + AR{zY + l, 
from which we can deduce the first terms of R: 

R{z) = z + ^z' + f-/ + 0{z') = zi^l + ^z' + + 0{z' 
Squaring R yields 
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Taking the reciprocal of the right-hand series finally yields 



p{z) = ^ (l - - + = ^ - - + 0{z% 

as requested. Thus, our fast algorithm to compute the coefficients C\ , . . . , Cfi IS clS follows: 

1. Compute R{z) mod ^2"+^ using the algorithm of g2ilwith G = Bt^ + Af^ + 1; 

2. Compute Q{z) = R{zf mod ^^"+5; 

3. Compute p{z) = l/Qiz) mod z^''+\ 

In the first step, we remark that our assumption R'{0) 7^ is indeed satisfied, hence 
R{z) mod z^"^*^ can be computed in 0(M(n)) operations, assuming 2,...,2n + 3 are 
units in K. Using the algorithm of ^2.11 the squaring and reciprocal necessary to recover 
p{z) mod admit the same complexity bound. This proves Theorem ^ 



4 Fast computation of isogenics 

In this section, we recall the basic properties of isogenics and an algorithm due to 
Elkies [221 that computes an isogeny of degree i in quadratic complexity 0(£^). Then, 
we design two fast variants of Elkies' algorithm, by exploiting the differential equations 
satisfied by some functions related to the Weierstrass function, proving Theorem |21 
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4.1 Isogenies 

The following properties are classical; all the ones not proved here can be found for 
instance in |49t i50j . Let E and E be two elliptic curves defined over K. An isogeny 
between E and £^ is a regular map I : E —>■ E that is also a group morphism. Hence, we 
have E ^ E/F, where F is the kernel of /; here, our isogenies are all non-zero. 

The most elementary example of an isogeny is the "multiplication by m" map which 
sends P & E to [m]P, where, as usual, the group law on E is written additively. If E is 
given through a Weierstrass model, the group law yields the following formulas for [m]P 
in terms of the Weber polynomials 4'm{x,y) [IHl P- 105]: 

Using the Weierstrass equation of E, the polynomial ipm{x, y) rewrites in terms of the so- 
called division polynomial fmix), which is univariate of degree 9(m^): ipm{x^y) = fm{x) 
if m is odd, ipm{x,y) = 2yfm{x) otherwise. 

Given an isogeny I : E E, there exist a unique isogeny (the dual isogeny) I : E ^ E 
and a unique integer i such that I o I = [^] ; the integer £ is called the degree of J; if I is 
separable, it equals the cardinality of its kernel. For instance, the degree of the isogeny [m] 
is and this is reflected by the degree of the division polynomials. 

Let E and E be two isogeneous elliptic curves in Weierstrass form, defined over K. 
Then the isogeny / between E and E can be written as 

I{x,y) = {Ux),cyr,{x)), (10) 

for some c in K. We say that / is normalized if the constant c equals 1; in this case, / 
is separable. We use an explicit form for such isogenies, extending results of Kohel 
§2.4] and Dewaghe |^ to the case of arbitrary degree i. 

Proposition 4.1. Let I : E ^ E be a normalized isogeny of degree i and let F be its 
kernel. Then I can be written as 

where D is the polynomial 



D{xy'' \D{x) 



D{x) = Y[{x-xq)= x'-^ - ax^-^ + (y2x'-^ - a^x'-^ + ■ ■ ■ (12) 

and N{x) is related to D{x) through the formula 

N(x) „ , , ,M(x) , , , , fD'{x)\' , , 

^ ' - ix-o- (Sx^ + A)-^ - 2{x? + Ax + S) — . (13) 



D{x) ' ' D[x) ' ' \ D{x 

Proof. Note first that given a subgroup F of -E(K), there can exist only one pair [E^ I) 
where E is in Weierstrass form and J is a normalized isogeny E ^ E having F as kernel. 
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In jni], Velu constructs the curve E and the normahzed isogeny /, starting from the 
coordinates of the points in its kernel F. A point P of coordinates {xp,yp) is sent by the 
isogeny / to a point of coordinates 

Xl(P) = Xp+ {Xp+Q - Xq) and ?//(p) = + ^ (yp+g - yg). 

Q£F* Q&F* 

From there, Velu uses the group law to get explicit expressions of the coordinates. More 
precisely, write F2 for the set of points in F that are of order 2. Then F can be written 

as 

F = {Oe} U F2 U Fodd U (-Fodd), 

where Fodd H (— Fodd) = and — Fodd denotes the set of opposite points of Fodd, so that 
D{x) rewrites as 

D{x) = YI{x-xq) Yl {x-XQf. 

Finally, let F~^ = F2 U Fodd- Then Velu gave the following explicit form for I{x,y) = 
(4(a;),y/;(x)): 

T ( ^ ^ V- f tQ , x^ + Axq + B 

h{x) = x+ }^ r— - + 4 



QeF+ 



X- Xq {x - XqY 



where Iq = 3xq + A if Q & F2 and tg = 2{3xq + A) otherwise. Observing that for 
Q G F2, Xq + Axq + B equals 0, one sees that Ix admits D for denominator, as claimed 
in Equation (fTTj) . 

Next, we split the sum over F+ into that for Q E F2 and that for Q G Fodd- The 
former rewrites as 

J2 + ^ 2^^! + AXQ + B 



^3 



X- Xq {x - Xq) 



2 



since these points satisfy Xq + Axq + B = 0; the sum for Q G Fodd rewrites as 



2 ^ \X — Xo (x — Xn) 



Therefore, we obtain 



, ^ f ^xl + A xl + AxQ + B \ 

h{x) = a; + > + 2 — — , 

q7^. V a; - Xq [x- Xq)2 J 

which can be rewritten as 

3x'^ + A x^ + Ax + B 
h[x) = x+ } [x-xq h2— — 

V X-Xq [x- Xq)^ 

This yields Equation (fT!?|) . □ 
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Though this is not required in what follows, let us mention how Vein's formulae enable 
one to construct the curve E. Let cr, a2, crs be as in Equation (|T^ and 

t = A{i-l) + 3{a^-2a2), 

w = 3Aa + 2B{i - 1) + 5(a^ - Saffa + Sag). 

Then the isogenous curve E has the Weiestrass equation = + AX + B, where 
A = A - 5t and B = B - 7w. 

The constant a introduced in the previous proposition is the sum of the abscissas of 
the points in the kernel F of /. In the important case where i is odd, the non-zero points 
in F come into pairs {{xq, i/q), {xq, —i/q)}, so that we will write 

D{x) = g{xf with g{x) = x^''^^'^ - qix^'~^^'^ + • ■ • , 

and a = 2qi. Then, we can replace D'{x)/ D{x) by 2g'{x)/g{x) in Proposition 14.11 



4.2 Elkies' quadratic algorithm 

From now on, we are given the two curves E and E through their Weierstrass equations, 
admitting a normalized isogeny I : E E oi degree ^. We will write 

E : y"^ = + Ax + B and E : y"^ = x^ + Ax + B. 

From this input, and possibly that of a, we want to determine the isogeny /, which we 
write as in Equation (jllj) 

(N{x) fN{x) 



We first describe an algorithm due to Elkies j22j, that we call Elkiesl998, whose com- 
plexity is quadratic in the degree i. In the next subsection, we give two fast variants of 
algorithm Elkiesl998, called fastElkies and fastElkies', of respective complexities 0{M{i)) 
and 0{M{e) logi). 

The algorithm Elkiesl998 was introduced for the prime degree case in [221, works 
for any i large enough. The first part of the algorithm aims at computing the expansion 
of N{x)/ D{x) at infinity; the second part amounts to recovering the power sums of the 
roots of D{x) from this expansion. 

To present these ideas, our starting remark is that the rational function N{x)/D{x) 
satisfies the non-linear differential equation 



(x' + Ax + B){ — f4 = — f4 ] +A{ — f4 + B. (14) 
^ \d{x)J \D{x)J \D{x)J ^ ^ 

This follows from Proposition 14.11 and the fact that / maps E onto E. Differentiating 
Equation p4|l leads to the following second-order equation: 

(3.^ , (E^) ' , ,A.,B) m) " - 3 f . i. ,15, 



D x) / ' ' \D x) / \D X 
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Writing the expansion of the rational function N{x)/D{x) at infinity 



N{x) , ^ hi 
— — — = x + } - 
D(x) ^ X' 



and identifying coefficients of x * from both sides of Equation yields the recurrence 



3 2k — 3 2(k — 3) 

hk = -r, n—, r > hjhk-T-i ; Ahun ; Bhk-^, for all k > 3, (16) 

{k-2){2k + 3) ^ ' ^ ' 2k + 3 2k + 3 - > v ; 



i=l 

with initial conditions 



hi = and /i2 ^ 



5 7 

The recurrence fll6|) is the basis of algorithm Elkiesl998; using it, one can compute 
/i3, . . . , hi_2 using 0(£^) operations in K. 

Elkies' algorithm Elkiesl998 assumes that a is given. Extracting coefficients in Equa- 
tion (|T!^ then yields 

hi = {2i + l)pi+i + {2i - l)Api^i + {2i - 2)Bpi^2, for all i > 1. (17) 

Since hi, ... , /i£_2 are known, P2, ■ ■ ■ ,P£-i can be deduced from the previous recurrence 
using 0{i) operations. The polynomial D{x) is then recovered, either by a quadratic 
algorithm or the faster algorithm of H2.2\ and N{x) is deduced using formula ()13|) . in 
0(M(£)) operations. 

This algorithm requires that 2, . . . , 2i—l be units in K. Its complexity is in 0(£^), the 
bottleneck being the computation of the coefficients hi, ... , hi-2- Observe the parallel 
with the computations presented in the previous section, where differentiating Weier- 
strass' equation yields the recurrence ((Tj), which appears as a particular case of the re- 
currence fll6|) (the former is obtained by taking A = B = in the latter). 



4.3 Fast algorithms 

We improve on the computation of the coefficients hi in algorithm Elkiesl998, the remain- 
ing part being unchanged. Unfortunately, we cannot directly apply the algorithm of 
to compute the expansion of N{x)/D{x) at infinity using the differential equation ()14|) 
since the equation obtained by the change of variables x ^ 1/x is singular at the origin 
To avoid this technical complication, we rather consider the power series 

S{x) = x + —j^^^ + ~Y4~^^ + 0{x^) e x + x^Kllx"^]] 



such that 



N{x) _ 1 
D(x) ^ 



2 ' 



remark that S satisfies the relation R = S o R, with the notation R{z) = 1/ a/ p{z) and 
R{z) = 1/ a/ p{z) introduced in 
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Applying the chain rule gives the following first order differential equation satisfied 
by S{x): 

{Bx^ + Ax^ + 1) S'{xf = l + A S{xY + B S{xf. 

Using this differential equation, we propose two algorithms to compute N{x)/D{x), de- 
pending on whether the coefficient a is known or not. In the algorithms, we write 

1 N(x) 
S{x) = xT(x^) and U{x) = G 1 + x^K[[x]] so that = xU 

The first algorithm, called fastElkies, assumes that a is known and goes as follows. 

1. Compute C{x) = {Bx^ + Ax^ + 1)-^ mod x^^-^ G K[[x]]; 

2. Compute S{x) mod x"^^ using the algorithm of ^2.41 with G{x,t) = C{x){l + At^ + 
Bt^), and deduce T(x) mod x^; 

3. Compute U{x) = 1/T(x)^ mod x^ using the algorithm in ^2.11 

4. Compute the coefficients hi, ... , /i£_2 of N{x)/D{x), using N{x)/ D{x) = xU (l/x); 

5. Compute the power sums P2, ■ ■ ■ ,Pe-i of -D(x), using the linear recurrence (|T7jl : 

6. Recover D{x) from its power sums, as described in ^2.2| 

7. Deduce N{x) using Equation (fT^ . 

Steps (1) and (5) have cost 0{i). Steps (2), (3), (6) and (7) can be performed in 0(M(£)) 
operations, and Step (4) requires no operation. This proves the first part of Theorem |^ 

For our second algorithm, that we call fastElkies', we do not assume prior knowl- 
edge of a. Its steps (l')-(3') are just a slight variation of Steps (l)-(3), of the same 
complexity 0(M(£)), up to constant factors. 

(r) Compute C{x) = {Bx^ + Ax^ + 1)-^ mod x^^-^ G K[[x]]; 

(2') Compute S{x) mod using the algorithm of gQlwith G{x,t) = C{x){l + At^ + 
Bt^), and deduce T(x) mod x^^"^; 

(3') Compute U{x) = 1/T(x)^ mod a;^^~^, using the algorithm in §2.1( 

(4') Reconstruct the rational function U{x); 

(5') Return N{x)/D{x) = xU{l/x). 

Using fast rational reconstruction. Step (4') can be performed in 0{M{i) logi) operations 
in K. Finally, it is easy to check that our algorithm fastElkies requires that 2, . . . , 2£ — 1 
be units in K, while algorithm fastElkies' requires that 2, . . . , 8£ — 5 be units in K. This 
completes the proof of Theorem |2l 

In the case of odd i, we can compute g{x) instead of D{x). Accordingly, we modify 
the recurrence relations, and compute fewer terms. Let qi,q2, . . . denote the power sums 
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of g{x), so that qi = Pi/2. Then, the coefficients hi and the power sums qi are related by 
the relation 

hi = {Ai + 2)gi+i + (4z - 2)Agi„i + (4i - A)Bqi_2. (18) 

To compute g{x) using algorithm fastElkies, it suffices to compute S{x) mod x^^^; then 
T{x) and U{x) are computed modulo x^^^^-*/^. Similarly, in algorithm fastElkies', it is 
enough to compute S{x) mod x^^, and T{x) and U{x) modulo x^^. 

5 Examples of isogeny computations 
5.1 Worked example 

Since the case of i odd is quite important in practice, we ffist give an example of such a 
situation (see below for an example with i = 6). Let 

E : y"^ = x^ + X + 1 and E : y"^ = x^ + 75x + 16 

be defined over Fioi, with i = 11 and a = 50. Since £ is odd, we will compute the 
polynomial g{x), which has degree 5. First, from the differential equation 

{x^ + x^ + l)S'{xy = 1 + 75S{xy + 16S{xf, S{0) = 0, ^'(0) = 1 

we infer the equalities 





C 


= 1 + 100 a; 


^ + 100 x*^ + x^ 


+ 2xi° + 0(x"). 




s 


= X + 68x^ 


+ 66 x^ + 60 x^ 


+ 84x1^ + 0(xi2 


so that 


T 


= 1 + 68x2 


+ 66 x^ + 60 x^ 


+ 84x5 + 0(x6). 


and 


T2 


= 1 + 35x2 


+ 31x3 + 98x^ 


+ 54x5 + 0(x''). 


whence 


U 


= 1 + 66x2 


+ 70x3 ^ iQ^4 


+ 96x^ + 0(x^). 



We deduce 



N(x) 66 70 16 96 ^ / 1 



At this stage, we know hi = 66, h2 = 70, h^ = 16, /i4 = 96, as well as qi = a/2 = 25. 
Equation (|T^ then writes 

hi - {At - 2)g,„i - {At - 4)g,_2 n i ^ ■ ^ ^ 

Qi+i = T--^ , for all 1 < z < 4 

At + 2 

and gives q2 = 43, q^ = 91, q^ = 86, q^ = 63. The main equation in ^2.21 writes 

5 ( (nr 43 2 91 3 86 4 63 5 

^ V X / " ^^^^ ( - ( 25 X + —x^ + y x'^ + — X* + Y x^ 

= exp6(76x + 29x2 + 37x2 + 29x^ + 48x^) , 

yielding g{x) = x^ + 76x^ + 89x3 _j_ 2/lx'^ + 97x + 5. For the sake of completeness, we have: 

N{x) = x" + 51x^° + 61x^ + 44x^ + 71x^ + 39x'' + 81x^ + 43x^ + 15x3 + 5x2 ^ 24x + 15. 

Had we computed the solution S'(x) at precision 0{x'^'^), the expansion at infinity of 
N{x) / D{x) would have been known at precision 0(l/x2^), and this would have sufficed 
to recover both A^(x) and D{x) by rational function reconstruction, without the prior 
knowledge of a. 
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5.2 Further examples. 

As it turns out, all the theory developed for the prime case in the SEA algorithm works 
also, mutatis mutandis, in the more general case of a cyclic isogeny of non prime degree. 
Consider the curve 

E -.y"^ = x^ + X + 3 

defined over Fioog- For £ = 6, we find that the modular poljTiomial of degree 6 (obtained 
as the resultant of the modular polynomials of degrees 2 and 3) has three roots, one of 
which is J = 248. Using the formulas of that are still valid, we find the isogenous 
curve 

E -.y^ = x^ + 830x + 82 
and a = 739, from which we obtain 

N{x) _x^ + 270 x^ + 325 x^ + 566 x^ + 382 x^ + 555 x + 203 
D{x) ~ x^ + 270 x^ + 289 x^ + 659 x^ + 533 x + 399 ' 

The denominator factors as 

(a;-66) (x - 23)^ (x - 818)^ 

The value x = 66 corresponds to one of the roots of x^ + x + 3 and is therefore the abscissa 
of a point of 2-torsion; 23 is the abscissa of a point of 3-torsion; 818 is the abscissa of a 
primitive point of 6-torsion. 

As an aside, let us illustrate the case of a non cyclic isogeny. The curve E happens 
to have rational 2-torsion; the subgroup E[2] of 2-torsion is non cyclic, being isomorphic 
to Z/2Z X Z/2Z. The denominator D{x) appearing in the isogeny I : E E = E/E[2] 
is simply D{x) = x^ + x + 3 and Equation (fT!?jl yields A^(x) = x^ + 1007x^ + 985 x + 1. 
From this, we can compute the equation of E, namely, y"^ = x^ + 16x + 192. 

6 A survey of previous algorithms for isogenies 

In this section, we recall and give complexity results for other known algorithms for 
computing isogenies. In what follows, we write the Weierstrass functions p and p of our 
two curves E and E as 

p{z) = ^ + X] ^ ~2 + 

i>l i>l 

All algorithms below require the knowledge of the expansion of these functions at least 
to precision i, so they only work under a hypothesis of the type p ^ i or p = 0. 

We can freely assume that these expansions are known. Indeed, by Theorem^ given 
A, B and A, B, we can precompute the coefficients q and q up to (typically) i = i — 1 
using 0(M(£)) operations in K, provided that the characteristic p of K is either or 
> £. This turns out to be negligible compared to the other costs involved in the following 
algorithms. 
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6.1 First algorithms 

A brute force approach to compute N{x)/ D{x) is to use the equation 

and the method of undetermined coefficients. This reduces to computing p{zy mod z^^~'^ 
for 1 < i < i and solving a hnear system with 2i — 1 unknowns. This direct method 
requires that 2, . . . , 4£ be units in K and its complexity is 0{i'^) operations in K, where 
2 < u; < 3 is the exponent of matrix multiplication. 

Another possible idea is to consider the rational functions N{x) / D{x) and N{x) / D{x) 
respectively associated to / and its dual /, noticing that by definition, 

N _ 

using the notation of Equation 0. However, algorithms for directly decomposing (j^e/i/jj ISHl 
123 Q] lead to an expensive solution in our case, since they require factoring the division 
polynomial fi, of degree 6(£^). Indeed, even using the best (sub- quadratic) algorithms 
for polynomial factorization [H2j, of exponent 1.815, this would yield an algorithm for 
computing isogenies of degree i in complexity more than cubic with respect to i, which 
is unacceptable. 

6.2 Stark's method 

To the best of our knowledge, the first subcubic method for finding A^ and D is due to 
Stark and amounts to expanding p as a continued fraction in p, using Equation (fTTHl . 
The fraction N/D is approximated by Pn/ln and the algorithm stops when the degree of 
g„ is £ — 1, yielding D. In particular, it works for any degree isogeny. Since p and p are 
in 1/ + K[[z^]], it is sufficient to work with series 'm Z = z^. 

1. T:= p{Z) + 0{Z% 

2. n := 1; 

3. go := 1; 

4. qi := 0; 

5. while deg(g„) < £ — 1 do 

{at this point, T{Z) = t.^Z""^ + ■ ■ ■ + + tiZ + ■ ■ ■ + 0(Z(^"d<=s9"-'^)-i)} 

(a) n := n + 1; 

(b) an ■■= 0; 

(c) while r > 1 do 

QjYi . (In ~t~ t—^Z , 



15 



T ■=T -t.rp'' = t-sZ-' + ■■■ ; 
r := s 

(d) g„ := a„g„_i + g„_2; 

(e) T := 1/T; 

6. Return D := qn- 

This algorithm (that we call Starkl972) requires 0{i) passes through Step (5); this 
bound is reached in general, with r = 1 at each step. The step that dominates the 
complexity is the computation of reciprocals in Step (5.e), with precision 2£— 1 — 2 degg„ — 
2r. The sum of these operations thus costs 0(£M(£)). The multiplications in Step (5.d) 
can be done in time 0{iM{i)) as well (these multiplications could be done faster if needed). 
Since the largest degree of the polynomials a„ is bounded by ^ — 1, computing all powers 
of p at Step (5.c) also fits within the 0(£M(£)) bound. Finally, knowing D{x), the 
numerator N{x) can be recovered in cost 0{M{i)) using Equation ()13|) . 

In the case where i is odd and we need g{x), as in the context of the SEA algorithm, 
we can compute it in 0(M(£)) operations by computing exp((logD)/2). 

To summarize, the total cost of algorithm Starkl972 is in 0{£M{i)). Remark that 
compared to the methods presented below, algorithm Starkl972 does not require the 
knowledge of a. Remark also that, even though r will be 1 in general, the computation 
of the powers p"^ in Step (5.c) could be amortized in the context of the SEA algorithm. 



6.3 Elkies' 1992 method 

We reproduce the method given in that we call Elkiesl992 (see also e.g., [TT| IHH]). 
We suppose that £ is odd, so that D{x) = g^x^, though minor modifications below would 
lead to the general solution. 

Differentiating twice Equation ^ yields 

= 120p3 + 72Ap + 485. 
More generally, we obtain equalities of the form 



d''p{z) 



fc + 1 



for some constants that satisfy the recurrence relation 

fik+i,j = (2j - 2)(2j - + (2j + l)(2j + 2)Afx,j+, + (2j + 2)(2j + 4)S/ife,,+2, 

with fik,k+i = (2 A; + 1)!. Using this recurrence relation, the coefficients iik,j, for k < d — 1 
and j < k + 1, can be computed in 0(£^) operations in K. 

Elkies then showed how to use these coefficients to recover the power sums q2, ■ ■ ■ ,qd 
of g, through the following equalities, holding for k > 1: 

{2k)\{ck - Cfc) = 2(/ifc,ogo H \- fik,k+iqk+i)- 
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Using these equalities, assuming that qi = a/2 and the coefficients Cfc, Ck and ^k,j are 
known, we can recover g2, • • • , Q'd by solving a triangular system, in complexity 0(£^). We 
can then recover g using either a quadratic algorithm, or the faster algorithm of §2.21 

There remains here the question whether the triangular system giving q2, . . . ,qd can 
be solved in quasi-linear time. To do so, one should exploit the structure of the triangular 
system, so as to avoid computing the 9(£^) constants fik,j explicitly. 

6.4 Atkin's method 

In 12], Atkin gave a formula enabling the computation of D{x) (see also |4()[ Formula 6.13] 
and (ISj) in the case where i is odd. We extend it, so as to cover the case of arbitrary i, 
this time recovering D[x). The equation we use is 

D{p{z)) = z'-^'eMm). (20) 

where 

Since d and the coefficients c^, Ck are all assumed to be known, one can deduce the series 
F{z) mod z^, provided that a is given. A direct method to determine D{x) is then to 
compute the exponential of F{z), and to recover the coefficients of D{x) one at a time, 
as shown in the following algorithm, called Atkinl992. As before, we use series in Z = . 

1. Compute the series Pi{Z) = p{Zy at order i, for 1 < i < i — 1] 

2. Compute G{Z) = exp^(F(Z)); 

3. T := G; 

4. D := 0; 

5. for i := i — 1 downto do 

{at this point, T = tZ"^ + ■ ■ ■ } 

(a) D := D + tz'] 

(b) T := T - tP,. 

Step (1) uses 0(£M(£)) operations; the cost of Step (2) is negligible, using either classical 
or fast exponentiation. Then, each pass through Step (5) costs 0{i) more operations, for 
a total of 0(£^). Thus, the total cost of this algorithm is in 0{iM{i)). If this algorithm 
is used in the context of the SEA algorithm. Step (1) can be amortized, since it depends 
on E only. Therefore, all the powers of p should be computed for the maximal value of 
i to be used, and stored. Hence, the cost of this algorithm would be dominated by that 
of Step (5), yielding a method of complexity 0(£^). 

A better algorithm for computing D{x), avoiding the computation of all powers of 
p{z), is based on the remark that Equation ()20j) rewrites 

d(-] =j2-2^((expoF)oJ), (21) 
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with T{x) = p ^(1/a;), where p ^ is the functional inverse of p. The expansion of 
at order 9(£) can be computed in 0{i) operations using the differential equation 

^'^""^'^ 4x{l + Ax^ + Bx^) °' ^'^''^ = ViTT + lFTW ^^^^ 
A linear differential equation follows: 

T{x) l + 3Ax^ + 4:Bx^ 



I{x) 2x{l + Ax^ + Bx^)' 

From there follows a linear differential equation for J7'(x) = x^^/^X(x) = X]i>o'^«'^*- 
Extracting coefficients in this equation then gives a linear recurrence 

a^+l = - 2(^+^lX2i + 3) ~ 3)5ai_2 + 2Ami„i) for i > 2, (23) 

with initial conditions oq = l,ai = 0,02 = This yields the following algorithm, 
called AtkinModComp: 

1. Compute G{Z) = exp^(F(Z)); 

2. Compute X(x) using Equation (j^ : 

3. Compute G{I) by modular composition (which is possible since G is in K[[Z]] = 
K[[x2]]); 

4. Deduce D using Equation (j2H). 

The cost of the algorithm is dominated by the composition of the series G = exp oF and 
I. From 3231 this can be done in 0{M{i)y/i + i'^) or 0{M{i) y/i \og£) operations in K. 

To do even better, it is fruitful to reconsider the series G{I) = (exp o F) oX used 
above, but rewriting it as exp o [F oX) instead; this change of point of view reveals close 
connections with our fastElkies algorithm. More precisely, Atkin's Equation ()20|) can be 
rewritten as 

D{p{x)) = exp (-ax^ + 2 ip{x) - p{x) 



We can then obtain D{l/x) as the following exponential: 

D (^] = exp (^-aJ2 + 2 j T j T (^-~{po J)(x)^ ^ (24) 



X 



Then, working out the details, the sequence of operations necessary to evaluate this 
exponential turns out to be the same as the one used in our algorithm fastElkies of §4.31 
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This does not come as a surprise: the relation ()17|1 used in our algorithm follows from 
formula |T!?|l . which can be rewritten as 

Ix-a- 2Vx3 + Ax + b(^x^ + Ax + B ^'^^^^ 



D{x) V D{x 

Then, Equation is nothing but an integral reformulation of this last equation, taking 
into account the fact that X satisfies the differential equation 



6.5 Summary 

In Tabled we gather the various algorithms discussed in this article, and compare these 
algorithms from two points of view: their complexity (expressed in number of operations 
in the base field K) and their need for a as input. 



algorithm 


complexity 


need of cr 


linear algebra 


o(r) 


no 


Starkl972 
Atkinl992 


o{m{i)) 

0(M(£)v^ + £'^) or 0(M(£)v/^log£) 

o{e) 
o{e) 

0(M(£)) 
0(M(£) log£) 


no 

yes 


AtkinModComp 
Elkiesl992 


yes 
yes 


Elkiesl998 
fastElkies 
fastElkies' 


yes 
yes 
no 



Table 1: Comparison of the algorithms 



7 Implementation and benchmarks 

We implemented our algorithms using the NTL C++ library j^lTTj and ran the program 
on an AMD 64 Processor 3400+ (2.4GHz). We begin with timings for computing the 
expansion of p, obtained over the finite field F102004+4683; they are given in Figure [H 

The shape of both curves indicates that the theoretical complexities - quadratic vs. 
nearly linear - are well respected in our implementation (note that the abrupt jumps 
at powers of 2 reflect the performance of NTL's FFT implementation of polynomial 
arithmetic). Moreover, the threshold beyond which our algorithm becomes useful over 
the quadratic one is reasonably small, making it interesting in practice very early. 

We now turn our attention to the pure isogeny part, concentrating on the case where 
£ is prime, in the context of the SEA algorithm. Hence, in this case, it suffices to compute 
the polynomial g{x) such that D{x) = g{xY. All algorithms can be adapted to make 
advantage of this simplification, as exemplified in ^4.31 for our algorithms fastElkies and 
fastElkies'. 

The first series of timings concerns the computation of isogenics over a small field, 
K = F1019+51, for the curve E : y"^ = x^ + 4589x + 91128. We compare in Figure El the 
performances of the algorithms Elkiesl992 from ^6. 31 and Elkiesl998 from ^4. 21 for isogenics 
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Figure 1: Timings for computing p on E : y = x + 4589x + 91128 over Fio2oo4_,_4683 

of moderate degree £ < 400. Figure El compares the timings obtained with the algorithm 
Elkiesl998 and our fast version fastElkies from ^4.H^ for isogenies of degree up to 6000. 




50 100 150 200 250 300 350 400 



Figure 2: Elkiesl992 vs. Elkiesl998. 





"isog.elkies98" 




"isog.laslelkies" 4/ 











Figure 3: Elkiesl998 vs. fastElkies. 



Next, we compare in Figure|3]the timings obtained by the 0{M{£)) algorithm fastElkies, 
that requires the knowledge of a, to those obtained by its 0{M{i)\ogi) counterpart 
fastElkies', that does not require this information. 



20 




2000 4000 BOOO 8D00 1D000 12D00 

Figure 4: FastElkies vs. FastElkies' 



In all figures, the degrees ^ of the isogenies are represented on the horizontal axis and 
the timings are given (in seconds) on the vertical axis. Again, the shape of both curves in 
Figure El shows that the theoretical complexities are well respected in our implementation. 
The curves in Figure|3]show that the theoretical ratio of log i between algorithms fastElkies 
and fastElkies' has a consequent practical impact. 

Next, in Tables |21 to |Hl we give detailed timings on computing ^-isogenies for the curve 

E : y"^ = x"^ + Ax + B 

where 

A = [10^^^°7rJ = 31415926 . . . 58133904, 

B = [10^^^°eJ = 27182818 . . . 94787610, 

for a few values of i, over the larger finite field Fio2oo4_,_4683, and using various meth- 
ods: algorithms Elkiesl992, Elkiesl998 and our fast variant fastElkies, Stark's algorithm 
Starkl972 and the two versions Atkinl992 and AtkinModComp of Atkin's algorithm. 

Tables El and El give timings for basic subroutines shared by some or all of the algo- 
rithms discussed. Table Elgives the timings necessary to compute the expansions of p and 
p, using either the classical algorithm or our faster variant: this is used in all algorithms, 
except our fastElkies algorithm. Table El gives timings for recovering g from its power 
sums, first using the classical quadratic algorithm, and then using fast exponentiation as 
described in ^2.2[ This is used in algorithms Elkiesl992, and Elkiesl998 and its variants. 

Tables 0] and El give the timings for algorithms Elkiesl992 on the one hand and 
Elkiesl998 and our variation fastElkies on the other hand. In Table I5 the columns /i 
and Pi give the time used to compute the coefficients fiij and the power sums pi. In 
Table El the column hi indicates the time used to compute the coefficients hi of the ra- 
tional function N/D, first using the original quadratic algorithm Elkiesl998, then using 
our faster variant fastElkies. The next column gives the time used to compute the power 
sums Pi from the hi using recurrence (|T8|l . 

Tables El and El give timings for our implementation of Atkin's original algorithm 
Atkinl992, as well as the faster version AtkinModComp using modular composition men- 
tioned in ^6.41 In Table El the column "exponential" compares the computation of exp(F) 
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Computing p and p 



order 



1013 
2039 
3019 
4001 
5021 



511 
1024 
1514 
2005 
2515 



quadratic 



8.6 
34.6 
75.7 
132.7 
209.3 



fast 



7.0 
29.9 
30.3 
31 
64.4 





Recovering g 


t 


quadratic 


last 


1013 




i.i 


2039 


17/1 

i (A 


l.D 


3019 






/inm 


66.9 


5.5 


5021 


106.2 


11.2 



Table 2: Computing p and p Table 3: Recovering g from its power sums 







Elkiesl992 




^ 


P 


/i 


Pi 


9 


1013 




10.4 


4.4 




2039 


See 


49.1 


17.9 


See 


3019 


Table H 


130.6 


38.9 


Table El 


4001 




263 


68.4 




5021 




496.5 


106.6 





Table 4: Algorithm Elkiesl992 





Elkiesl998 and fastElkies 


^ 


hi 




Pi 


9 




quadratic 


fast 






1013 


4.4 


4.5 


0.05 




2039 


17.3 


9.6 


0.1 


See 


3019 


38.0 


19.5 


0.16 


Table IHl 


4001 


67.2 


20.0 


0.21 




5021 


105.0 


40.7 


0.27 





Table 5: Algorithms Elkiesl998 and fastElkies 



using the naive exponentiation algorithm to the computation using the faster algorithm 
presented in ^2.2j the column p'^ gives the time for computing all the series p{zY and 
the column g that for recovering the coefficients of g from its power sums. Table [3 
gives timings obtained using the two modular composition algorithms mentioned in ^2.5t 
called here ModCompl and ModComp2; the previous columns give the time for computing 
exp(F) and that for computing the requested power of X; the last column gives the time 
to perform the final multiplication. 

Asymptotically, algorithm ModComp2 is faster than algorithm ModCompl, so that the 
timings in Table[7|might come as a surprise. The explanation is that, for the problem sizes 
we are interested in, the predominant step of algorithm ModCompl is the one based on 
polynomial operations, while the step based on linear algebra operations takes only about 
10% of the whole computing time. Thus, the practical complexity of this algorithm in the 
considered range (1000 < I < 6000) is proportional to M(£)v^, while that of algorithm 
ModComp2 is proportional to M{€)\/T\ogI. Moreover, the proportionality constant is 
smaller in the built-in NTL function performing ModCompl than in our implementation 
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Algorithm Aktinl992 




n 
1 




exponential 




9 






naive 


last 






1013 




O O A 

88.4 


"1 o 

1.2 


TO Q 

/ z.o 


A A 

4.4 


or\ on 

2039 


See 


o ( U. i 




oU4.y 


17.7 


on 1 n 

3019 




955.9 


5.1 


1 tJtJ.O 


o o n 

38.9 


/inni 

4UUi 




1503 


5.2 


1218.9 


D ( .0 


5021 




3180 


10.8 


2506.4 


108.7 



Table 6: Atkin's original algorithm, variations for exp(F) 







Algorithm AtkinModComp 




i 


P, P 


exp(F) 


ji-i 


modular composition 


9 










ModCompl 


ModComp2 




1013 




1.2 


2.7 


14.3 


35.6 


0.2 


2039 


See 


2.5 


6.6 


45.8 


111.9 


0.4 


3019 


Table 12 


5.1 


10.4 


95.3 


241 


0.7 


4001 




5.2 


11.6 


143.2 


338 


0.9 


5021 




10.9 


20.9 


240 


642 


1.4 



Table 7: Atkin's algorithm with modular composition 



of ModComp2. 

Notice that in all the columns labelled "fast" in Tables I2HZI the timings reflect the 
already mentioned (piecewisely almost constant) behaviour of the FFT: polynomial multi- 
plication in the degree range 1024-2047 is roughly twice as fast as in the range 2047-4095 
and roughly four times as fast as in the range 4096-8191. 

Finally, TablelHlgives timings for Stark's algorithm Starkl972; apart from the common 
computation of p and p, we distinguish the time necessary to compute all inverses (the 
quadratic algorithm when available, followed by that using fast inversion) and that for 
deducing the polynomials g„. 



i 


P, P 


Inverses 


Qn 






quadratic 


fast 




1013 




23542 


1222.7 


28.0 


2039 


See 


> 100000 


5113.4 


116.9 


3019 


Table El 




12182 


258 


4001 






20388 


418.6 


5021 






38910 


663.1 



Table 8: Stark's algorithm Starkl972 
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Conclusion 



The complexity analyses of the algorithms we have surveyed shows that for the case of 
a large prime characteristic and for a reasonably large degree C. of the isogeny, our new 
0(M(£)) algorithm improves over previously known techniques. 

The current implementation of our algorithm can be further optimized to make it the 
algorithm of choice for smaller values of the degree. Indeed, it is known that algorithms 
based on Newton iteration present certain redundancies (coefficients that can be predicted 
in advance, repeated multiplicands). Removing these redundancies is feasible (see 
allowing one to achieve constant-factor speed-ups. For the moment, our implementation 
relies only partially on these techniques; we believe that further programming effort would 
bring practical improvements by non-negligible constant factors. 

Another direction for future work is to adapt our ideas to the case of a small charac- 
teristic. In this respect, modifying the last phase of the algorithm of Joux and Lercier 
[31 seems a promising search path. 

Acknowledgments. We thank Pierrick Gaudry for his remarks during the elaboration 
of the ideas contained in this work. 
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